function t = targetfun(param,mw,mw16,p)

% replace parameters to be calibrated internally
p.m = param(1);
p.n = param(2);
p.a = param(3);
p.b = param(4);
p.ak = param(5);

if numel(param)>6 % scale parameter A is provided
    p.A = param(7);
end

% compute q from qs, i.e., param(6)
% first compute qinfinity:
p.q = 0; % some value needed to compute next line, but the exact value doesn't matter
qinfinity = log(integral(@(x) exp((logcapprod(x,p)-p.q)*(p.lambda-1)),0,1)^(1/(1-p.lambda)));
qmax = qinfinity + log(0.21)-2*p.qtol; % maximal q in 1980 s.t. q in 2016 is below q_infinity-p.qtol
% next qzero (w/o minimum wage):
rno = laborassignment(0,0,0,0,p); % equilibrium without capital
logomega = log(rno.w) + logcapprod(rno.x,p) - p.q...
    - loglaborprod(rno.sgrid,rno.x,p); % effective cost of labor relative to capital
qzero = -max(logomega);
% compute 1980 q
q = min(qmax,qzero + param(6)*(qmax-qzero));

% solve model for 1980 q and minimum wage
r1980 = solvemodelmw(p,q,mw);

% compute q in 2016 from q in 1980, using the indicated price index
% see "Acemoglu Loebbing Capital Shares and User Cost.xlsx" for price
% indices
if p.qindex == "eqqa" % non-residential equipment, quality-adjusted
    q16 = q - log(0.15);
elseif p.qindex == "eq" % non-residential equipment
    q16 = q - log(0.42);
elseif p.qindex == "sw" % software
    q16 = q - log(0.26);
else % baseline: non-residential equipment and software, quality-adjusted
    q16 = q - log(0.21);
end

% solve model for 2016 q and minimum wage
r2016 = solvemodelmw(p,q16,mw16);

% compute skill indices of 10, 30, 50, 90 percentiles of 1980 wage
% distribution
s10 = 0.1*(1-r1980.sbar)+r1980.sbar;
s30 = 0.3*(1-r1980.sbar)+r1980.sbar;
s50 = 0.5*(1-r1980.sbar)+r1980.sbar;
s90 = 0.9*(1-r1980.sbar)+r1980.sbar;

% compute calibration targets
if r1980.exitflag == 1
    % 10, 30 and 90 percentiles of 1980 wage distribution
    ind10 = find(r1980.sgrid>s10,1); % index of p10
    w10 = interp1([r1980.sgrid(ind10-1) r1980.sgrid(ind10)],[r1980.w(ind10-1) r1980.w(ind10)],s10);
    ind30 = find(r1980.sgrid>s30,1); % index of p30
    w30 = interp1([r1980.sgrid(ind30-1) r1980.sgrid(ind30)],[r1980.w(ind30-1) r1980.w(ind30)],s30);
    ind50 = find(r1980.sgrid>s50,1); % index of p50
    w50 = interp1([r1980.sgrid(ind50-1) r1980.sgrid(ind50)],[r1980.w(ind50-1) r1980.w(ind50)],s50);
    ind90 = find(r1980.sgrid>s90,1); % index of p90
    w90 = interp1([r1980.sgrid(ind90-1) r1980.sgrid(ind90)],[r1980.w(ind90-1) r1980.w(ind90)],s90);
    
    % log of 50 percentile of wages in 1980
    logw50 = log(w50);
    
    % log of 30-10 percentile ratio of wages in 1980
    logw3010 = log(w30)-log(w10);
    
    % log of 50-30 percentile ratio of wages in 1980
    logw5030 = log(w50)-log(w30);
    
    % log of 50-10 percentile ratio of wages in 1980
    logw5010 = log(w50)-log(w10);
    
    % log of 90-50 percentile ratio of wages in 1980
    logw9050 = log(w90)-log(w50);

    % factor income of equipment capital and software relative to labor in 1980
    alphak = r1980.alphak;
    
    % log average wage in 1980
    logw = log((1-alphak)*r1980.y);
    
else
   
    logw50 = NaN;
    logw3010 = NaN;
    logw5030 = NaN;
    logw5010 = NaN;
    logw9050 = NaN;
    alphak = NaN;
    logw = NaN;
    
end
   
if r2016.exitflag == 1
    % 10, 30 and 90 percentiles of hypothetical 2016 wage distribution
    ind10new = find(r2016.sgrid>s10,1); % index of p10
    w10new = interp1([r2016.sgrid(ind10new-1) r2016.sgrid(ind10new)],[r2016.w(ind10new-1) r2016.w(ind10new)],s10);
    ind30new = find(r2016.sgrid>s30,1); % index of p30
    w30new = interp1([r2016.sgrid(ind30new-1) r2016.sgrid(ind30new)],[r2016.w(ind30new-1) r2016.w(ind30new)],s30);
    ind50new = find(r2016.sgrid>s50,1); % index of p50
    w50new = interp1([r2016.sgrid(ind50new-1) r2016.sgrid(ind50new)],[r2016.w(ind50new-1) r2016.w(ind50new)],s50); 
    ind90new = find(r2016.sgrid>s90,1); % index of p90
    w90new = interp1([r2016.sgrid(ind90new-1) r2016.sgrid(ind90new)],[r2016.w(ind90new-1) r2016.w(ind90new)],s90);

    % change in 10 percentile log wage 1980-2016 due to automation
    deltalogw50 = log(w50new)-logw50;
    
    % change in log of 30-10 percentile ratio of wages 1980-2016 due to
    % automation
    deltalogw3010 = log(w30new) - log(w10new) - logw3010;
    
    % change in log of 50-30 percentile ratio of wages 1980-2016 due to
    % automation
    deltalogw5030 = log(w50new) - log(w30new) - logw5030;
    
    % change in log of 90-50 percentile ratio of wages 1980-2015 due to
    % automation
    deltalogw9050 = log(w90new) - log(w50new) - logw9050;
    
    % factor income of equipment and software relative to labor in
    % hypothetical 2016
    alphaknew = r2016.alphak;
    
    % change in factor income share of equipment and software 1980-2016 due
    % to automation
    deltaalphak = alphaknew - alphak;
    
    % average wage in hypothetical 2016
    logwnew = log((1-alphaknew)*r2016.y);
    
    % change in average wages 1980-2016 due to automation
    deltalogw = logwnew - logw;
    
else
    
    deltalogw50 = NaN;
    deltalogw3010 = NaN;
    deltalogw5030 = NaN;
    deltalogw9050 = NaN;
    deltaalphak = NaN;
    deltalogw = NaN;
    
end

% collect results

t.targets = [logw5010; logw9050; alphak; deltalogw3010; deltalogw5030; deltalogw9050; deltaalphak; logw50; deltalogw; deltalogw50];
t.efs = [r1980.exitflag; r2016.exitflag];
t.details1980 = r1980;
t.details2016 = r2016;

end